library(data.table)
library(Seurat)
library(harmony)
library(SingleR)
library(tidyseurat)
library(tidyverse)
source('00_util_scripts/mod_seurat.R')

# read raw data -------------
count_mtx <- fread('mission/FPP/psoriasis/GSE162183_Raw_gene_counts_matrix.tab.gz')

count_mtx[1:5, 1:5]

column_to_rownames(count_mtx, 'V1') -> count_mtx

sobj <- CreateSeuratObject(count_mtx, min.cells = 3, min.features = 200)

sobj <- sobj |>
  PercentageFeatureSet('^MT-', col.name = 'mito.ratio')

sobj |> VlnPlot('mito.ratio', pt.size = 0)

sobj <- sobj |>
  filter(mito.ratio < 10)

sobj <- sobj |>
  mutate(group = str_extract(orig.ident, 'Ctrl|Psor'))

# pre-process -------------
sobj <- sobj |>
  quick_process_seurat()

# identify cell type ---------
hpca <- celldex::HumanPrimaryCellAtlasData()

sobj <- sobj |>
  mark_cell_type_singler(hpca, new_label = 'hpca_main')

# save and read -----------
write_rds(sobj, 'mission/FPP/psoriasis/gao2021harbin.rds')

sobj |>
  filter(str_detect(hpca_main, 'Kera')) |>
  write_rds('mission/FPP/psoriasis/gao2021kera.rds')

sobj.kera <- read_rds('mission/FPP/psoriasis/gao2021kera.rds')

# visualization mission -----------
kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps')|>
  str_to_upper()

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

kegg_d.mva <-
  c('FDFT1','LSS','SQLE','GGPS1','DOLK','PDSS1','PDSS2') |>
  fct_inorder()

key_cytokine <- c('Il6','Ccl20','Tslp','Flt3l','Csf2','Tnf') |>
  str_to_upper()

celltyps.list <- sobj$hpca_main |> unique()

celltyps.list <- celltyps.list |>
  set_names(celltyps.list)

Idents(sobj) <- 'hpca_main'

celltyps.psor.vs.ctrl <- celltyps.list |>
  map(\(x)sobj |>
        FindMarkers(ident.1 = 'Psor', group.by = 'group', subset.ident = x) |>
        as_tibble(rownames = 'gene')) |>
  list_rbind(names_to = 'cell.type')

celltyps.psor.vs.ctrl <- celltyps.psor.vs.ctrl |>
  filter(p_val_adj < .05)

celltyps.psor.vs.ctrl |>
  right_join(mva_fct) |>
  filter(!is.na(avg_log2FC)) |>
  ggplot(aes(ordered, cell.type, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_color_distiller(palette = 'RdBu',limits = c(-4,4)) +
  theme_pubr(x.text.angle = 45) +
  labs(x = 'gene')

celltyps.psor.vs.ctrl |>
  filter(gene %in% key_cytokine) |>
  ggplot(aes(gene, cell.type, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_color_distiller(palette = 'RdBu',limits = c(-4,4)) +
  theme_pubr(x.text.angle = 45) +
  labs(x = 'gene')

sobj.kera |>
  get_abundance_sc_long('CCL20') |>
  mutate(group = ifelse(str_detect(.cell, 'Ctrl'), 'Ctrl', 'Psoriasis')) |>
  ggplot(aes(group, .abundance_RNA, fill = group)) +
  geom_violin() +
  stat_summary(geom = 'crossbar', fun = 'mean',
               color = 'red', width = .3) +
  theme_pubr(legend = 'none') +
  labs(title = 'CCL20 RNA expression in keratinocytes',
       y = 'Normalized expression')

sobj.kera <- sobj.kera |>
  quick_process_seurat(skip_norm = T)

sobj.kera |>
  FindAllMarkers(features = 'TRPV3', only.pos = T)

sobj.kera <- sobj.kera |>
  mutate(trpv3.hi = ifelse(seurat_clusters %in% c(2,4,6,8),
                           'TRPV3_high', 'TRPV3_low'))

sobj.kera |>
  VlnPlot(features = 'TRPV3', group.by = 'trpv3.hi', pt.size = 0) +
  stat_summary(geom = 'crossbar', fun = 'mean',
               color = 'red', width = .3) +
  theme_pubr(legend = 'none') +
  labs(title = 'CCL20 RNA expression in keratinocytes',
       y = 'Normalized expression', x = 'Subtype')

sobj.kera |>
  DimPlot(group.by = 'trpv3.hi') +
  ggtitle('Keratinocyte subtype on TRPV3 expression')

sobj.kera <- sobj.kera |>
  mutate(subgroup = str_c(trpv3.hi,'_', group)) 

sobj.kera |>
  DotPlot(kegg_mva, group.by = 'subgroup') +
  theme_pubr(x.text.angle = 45, legend = 'right') +
  labs(y = 'Subgroup', x = 'Gene',
       title = 'MVA pathway in TRPV3-hi/lo keratinocyte in psoriasis/ctrl skin')

sobj.kera |>
  DotPlot(key_cytokine, group.by = 'subgroup') +
  theme_pubr(x.text.angle = 45, legend = 'right') +
  labs(y = 'Subgroup', x = 'Gene',
       title = 'Cytokine in TRPV3-hi/lo keratinocyte in psoriasis/ctrl skin')
